Use R package UScensus2010county to complete the following tasks: (20 pt.)
Plot a map of New York counties using the plot function.
if (!require(UScensus2010)) {install.packages("UScensus2010")}
library(UScensus2010)
if (!require(UScensus2010county)) {install.county("osx")}
if (!require(UScensus2010tract)) {install.tract("osx")}
library(UScensus2010county)
library(UScensus2010tract)
data(new_york.county10)
shp <- new_york.county10
df <- shp@data
df
plot(shp)
Plot a map of New York counties using the qtm function.
library(UScensus2010county)
data(new_york.county10)
shp <- new_york.county10
df <- shp@data
df
plot(shp)
if (!require(tmap)) {install.packages("tmap")}
library(tmap)
qtm(shp)
How many counties in New York State?
shp <- new_york.county10
df <- shp@data
nrow(df)
[1] 62
What’s the 3-digit fips code of Broome County?
shp <- new_york.county10
df <- shp@data
sel <- df$fips[df$NAME10 == "Broome"]
sel
[1] "36007"
sprintf('%03d', as.numeric(sel) %% 1000)
[1] "007"
Compute descriptive statistics of the population column (P0010001), including total, minimum, maximum, mean, median, and skewness.
library(UScensus2010county)
data(new_york.county10)
shp <- new_york.county10
df <- shp@data
pop <- df$P0010001
sum(pop)
[1] 19378102
min(pop)
[1] 4836
max(pop)
[1] 2504700
mean(pop)
[1] 312550
median(pop)
[1] 91301
library(e1071)
skewness(pop)
[1] 2.517639
Create a histogram and a boxplot of the population.
hist(pop)
boxplot(pop)
Use R package UScensus2010tract to complete the following tasks: (20 pt.)
if (!require(UScensus2010tract)) {install.tract("osx")}
library(UScensus2010tract)
Plot a map of New York census tracts using the plot function.
data("new_york.tract10")
shp <- new_york.tract10
df <- shp@data
df
plot(shp)
Compute the total population based on census tracts.
data("new_york.tract10")
shp <- new_york.tract10
df <- shp@data
pop <- sum(df$P0010001)
pop
[1] 19378102
Select all census tracts in Broome County and plot the map.
data("new_york.tract10")
shp <- new_york.tract10
df <- shp@data
sel <- df$county == "007"
plot(shp[sel,])
What’s the total population of Broome County?
Broomepop <- df$P0010001[df$county == "007"]
sum(Broomepop)
[1] 200600
Create a histogram and a boxplot of population based on census tracts of Broome County.
hist(Broomepop)
boxplot(Broomepop)
Select the first five columns of the shapefile of Broome County census tract; add a new population ratio column (= census tract population / county population); save the new shapefile to the hard drive.
dfnew <- df[,1:5]
shp@data <- dfnew
shp@data
shp$ratio <- shp$P0010001/sum(Broomepop)
writeOGR(shp,dsn = "output",layer = "newshp",driver = "ESRI Shapefile",overwrite_layer = TRUE)
NAs introduced by coercion
newshp <- readOGR(dsn = "output", layer = "newshp")
OGR data source with driver: ESRI Shapefile
Source: "C:\Users\Yurui\Documents\GEOG_533\Lab09\GEOG-533\Labs\Lab_09\output", layer: "newshp"
with 4919 features
It has 6 fields
newshp@data
Use R packages raster and ncdf4 to complete the following tasks: (20 pt.)
if (!require(ncdf4)) {install.packages("ncdf4")}
library(ncdf4)
library(raster)
Load the multi-band raster dataset “NDVI.nc” into R.
ndvi <- brick("NDVI.nc")
ndvi
class : RasterBrick
dimensions : 360, 720, 259200, 36 (nrow, ncol, ncell, nlayers)
resolution : 0.5, 0.5 (x, y)
extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
data source : C:\Users\Yurui\Documents\GEOG_533\Lab09\GEOG-533\Labs\Lab_09\NDVI.nc
names : X2000.01.01, X2000.02.01, X2000.03.01, X2000.04.01, X2000.05.01, X2000.06.01, X2000.07.01, X2000.08.01, X2000.09.01, X2000.10.01, X2000.11.01, X2000.12.01, X2001.01.01, X2001.02.01, X2001.03.01, ...
Date : 2000-01-01, 2002-12-01 (min, max)
varname : NDVI
Get the basic information about the dataset, including the number of rows, columns, cells, and bands; spatial resolution, extent, bounding box, and projection.
nrow(ndvi)
[1] 360
ncol(ndvi)
[1] 720
ncell(ndvi)
[1] 259200
extent(ndvi)
class : Extent
xmin : -180
xmax : 180
ymin : -90
ymax : 90
bbox(ndvi)
min max
s1 -180 180
s2 -90 90
res(ndvi)
[1] 0.5 0.5
projection(ndvi)
[1] "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"
Aggregate all bands to generate a mean NDVI raster and a maximum NDVI raster; save the two new raster datasets to the hard drive.
ndvimean <- mean(ndvi)
ndvimean
class : RasterLayer
dimensions : 360, 720, 259200 (nrow, ncol, ncell)
resolution : 0.5, 0.5 (x, y)
extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
data source : in memory
names : layer
values : 0, 0.8537528 (min, max)
ndvimax <- max(ndvi)
ndvimax
class : RasterLayer
dimensions : 360, 720, 259200 (nrow, ncol, ncell)
resolution : 0.5, 0.5 (x, y)
extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
data source : in memory
names : layer
values : 0, 0.9681 (min, max)
writeRaster(ndvimean, "output/ndvimean.nc", overwrite = TRUE)
class : RasterLayer
dimensions : 360, 720, 259200 (nrow, ncol, ncell)
resolution : 0.5, 0.5 (x, y)
extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
data source : C:\Users\Yurui\Documents\GEOG_533\Lab09\GEOG-533\Labs\Lab_09\output\ndvimean.nc
names : layer
zvar : layer
writeRaster(ndvimax, "output/ndvimax.nc", overwrite = TRUE)
class : RasterLayer
dimensions : 360, 720, 259200 (nrow, ncol, ncell)
resolution : 0.5, 0.5 (x, y)
extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
data source : C:\Users\Yurui\Documents\GEOG_533\Lab09\GEOG-533\Labs\Lab_09\output\ndvimax.nc
names : layer
zvar : layer
Plot the maps of monthly NDVI of the year 2001
ndvi2001 <- ndvi[[13:24]]
plot(ndvi2001)
Create histograms of monthly NDVI of the year 2001.
hist(ndvi2001)
size changed to n because it cannot be larger than n when replace is FALSE39% of the raster cells were used (of which 42% were NA). 58093 values used.
Plot the NDVI map of July 2000; click any location with data on the map and retrieve the NDVI time series for all years; plot the NDVI time series of the selected location.
plot(ndvi, 7)
values <- click(ndvi, n=1, xy=FALSE)
plot(as.vector(values), type = "b")
no non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -InfError in plot.window(...) : need finite 'xlim' values